#include "TwoDimProblem.h"
using namespace std;
using namespace Utilities;

TwoDimProblem::TwoDimProblem(string initFile)
{
	this->paraFile = initFile;
	FileOperator paraFileOperator = FileOperator(paraFile);
	paraFileOperator.ReadTwoDimProblemInfo(gridNumX, gridNumY, leftBoundary, rightBoudary, upBoundary, downBoundary);
	normalizeBoundaryConditionValues();
	
	this->vertexNumX = gridNumX + 1;
	this->vertexNumY = gridNumY + 1;

	meshX.resize(gridNumX);
	meshY.resize(gridNumY);

	coordX.Resize(vertexNumX, vertexNumY);
	coordY.Resize(vertexNumX, vertexNumY);
	potential.Resize(vertexNumX, vertexNumY);
	elecFieldX.Resize(vertexNumX, vertexNumY);
	elecFieldY.Resize(vertexNumX, vertexNumY);
	elecFieldAbs.Resize(vertexNumX, vertexNumY);
	chargeAtVertex.Resize(vertexNumX, vertexNumY);
	
	rhs.resize(vertexNumX * vertexNumY);

	chargeConc.Resize(gridNumX, gridNumY);
	rDielectric.Resize(gridNumX, gridNumY);
	deltaX.Resize(gridNumX, gridNumY);
	deltaY.Resize(gridNumX, gridNumY);

	initAndNormPhysicalParameter();
}


TwoDimProblem::~TwoDimProblem(void)
{
}

void TwoDimProblem::Solve()
{
	solvePotential();
}

void TwoDimProblem::OutputResult()
{
	ParameterNormalization theNorm = ParameterNormalization();
	//theNorm.ConversePotentialVector(potential.DataVector, potential.DataVector, Pull);
	//PrintConsole::PrintMy2DVector(potential);
	vector<My2DVector<double>> output2DVectors;
	output2DVectors.push_back(coordX);
	output2DVectors.push_back(coordY);
	output2DVectors.push_back(potential);

	theNorm.ConverseLengthVector(output2DVectors[0].DataVector, output2DVectors[0].DataVector, Pull);
	theNorm.ConverseLengthVector(output2DVectors[1].DataVector, output2DVectors[1].DataVector, Pull);
	theNorm.ConversePotentialVector(output2DVectors[2].DataVector, output2DVectors[2].DataVector, Pull);

	FileOperator outFileOrigin = FileOperator("TwoDimOut_Origin.txt");
	string title = "X length\tY length\tPotential";
	outFileOrigin.Output2DVectorsForOrigin(title, output2DVectors);
 
	FileOperator outFileMatlab = FileOperator("TwoDimOut_Matlab.txt");
	title = "X coordinate";
	outFileMatlab.Output2DVectorForMatlab(title, output2DVectors[0], true);
	title = "Y coordinate";
	outFileMatlab.Output2DVectorForMatlab(title, output2DVectors[1], false);
	title = "Potential";
	outFileMatlab.Output2DVectorForMatlab(title, output2DVectors[2], false);
}

void TwoDimProblem::initAndNormPhysicalParameter()
{
	FileOperator paraFileOperator = FileOperator(this->paraFile);
	paraFileOperator.ReadTwoDimProblem(meshX, meshY, rDielectric, chargeConc);

	//for data monitoring
	//PrintConsole::PrintVector(meshX);
	//PrintConsole::PrintVector(meshY);
	
	ParameterNormalization theNormalization = ParameterNormalization();
	//the internal vector is used in the parameter converson
	theNormalization.ConverseLengthVector(meshX, meshX, Push);
	theNormalization.ConverseLengthVector(meshY, meshY, Push);
	theNormalization.ConveseConcVector(chargeConc.DataVector, chargeConc.DataVector, Push);

	//Fill in the (Nx+1)*(Ny+1) values
	//Fill in the coordX and coordY
	for (int iy = 0; iy != coordX.ColSize; ++iy )
	{
		for (int ix = 0; ix != coordX.RowSize; ++ix)
		{
			//bracket is important
			coordX(ix, iy) = ((ix == 0) ? 0 : coordX(ix-1, iy))
				           + ((ix == 0) ? 0 : meshX[ix-1]);
		}
	}
	for (int ix = 0; ix != coordY.RowSize; ++ix)
	{
		for (int iy = 0; iy != coordY.ColSize; ++iy)
		{
			coordY(ix, iy) = ((iy == 0) ? 0 : coordY(ix, iy-1))
				           + ((iy == 0) ? 0 : meshY[iy-1]);
		}
	}

	//for data monitoring
	//PrintConsole::PrintMy2DVector(coordX);
	//PrintConsole::PrintMy2DVector(coordY);

	//fill in charge at vertex
	for (int ix = 0; ix != chargeAtVertex.RowSize; ++ix)
	{
		for (int iy = 0; iy != chargeAtVertex.ColSize; ++iy)
		{
			chargeAtVertex(ix, iy) = ((ix-1>=0                     && iy<chargeAtVertex.ColSize-1) ? chargeConc(ix-1, iy)   * meshX[ix-1] * meshY[iy]   * 1.0/4 : 0)  //up-left
								   + ((ix<chargeAtVertex.RowSize-1 && iy<chargeAtVertex.ColSize-1) ? chargeConc(ix, iy)     * meshX[ix]   * meshY[iy]   * 1.0/4 : 0)  //up-right
								   + ((ix-1>=0                     && iy-1>=0)                     ? chargeConc(ix-1, iy-1) * meshX[ix-1] * meshY[iy-1] * 1.0/4 : 0)  //down-left
								   + ((ix<chargeAtVertex.RowSize-1 && iy-1>=0)                     ? chargeConc(ix, iy-1)   * meshX[ix]   * meshY[iy-1] * 1.0/4 : 0);  //down-right; 
		}
	}
	//data monitoring
	//PrintConsole::PrintMy2DVector(chargeAtVertex);
}

void TwoDimProblem::normalizeBoundaryConditionValues()
{
	ParameterNormalization theNormalization = ParameterNormalization();
	//in artificial boundary conditions, the value is 0, so there is no need to convere it.
	switch (leftBoundary.Type)
	{
	case First:
		leftBoundary.ValuePotential = theNormalization.PushPotential(leftBoundary.ValuePotential);
		break;
	case Second:
		leftBoundary.ValueElecFieldX = theNormalization.PushElecField(leftBoundary.ValueElecFieldX);
		break;
	}

	switch (rightBoudary.Type)
	{
	case First:
		rightBoudary.ValuePotential = theNormalization.PushPotential(rightBoudary.ValuePotential);
		break;
	case Second:
		rightBoudary.ValueElecFieldX = theNormalization.PushElecField(rightBoudary.ValueElecFieldX);
		break;
	}

	switch (upBoundary.Type)
	{
	case First:
		upBoundary.ValuePotential = theNormalization.PushPotential(upBoundary.ValuePotential);
		break;
	case Second:
		upBoundary.ValueElecFieldY = theNormalization.PushElecField(upBoundary.ValueElecFieldY);
		break;
	}

	switch (downBoundary.Type)
	{
	case First:
		downBoundary.ValuePotential = theNormalization.PushPotential(downBoundary.ValuePotential);
		break;
	case Second:
		downBoundary.ValueElecFieldY = theNormalization.PushElecField(downBoundary.ValueElecFieldY);
		break;
	}
}

void TwoDimProblem::initRHSVector()
{
	int index = 0;
	for (int ix = 0; ix != vertexNumX; ++ix)
	{
		for (int iy = 0; iy != vertexNumY; ++iy)
		{
			index = chargeAtVertex.IndexAt(ix, iy);
			rhs[index] = - chargeAtVertex(ix, iy);
		}
	}
	//for data monitoring
	//PrintConsole::PrintVector(rhs);
}

void TwoDimProblem::refreshRHSVector()
{
	int ix = 0;
	int iy = 0;
	int index = 0;
	My2DVector<double> positoner = My2DVector<double>(vertexNumX, vertexNumY);
	//for Artificial boundary conditions, the value is 0, so no change happens to the right-hand side
	//firstly judge if it is the second boundary condition and after that judge if it is the first boundary condition.
	//it is because first boundary condition is prior to second boundary conditions, so refresh it last to override related right-hand side
	
	//////////////////// for Neumann boundary conditions //////////////////
	if (leftBoundary.Type == Second)
	{
		ix = 0;
		for (iy = 0; iy != vertexNumY; ++iy)
		{
			index = positoner.IndexAt(ix, iy);
			rhs[index] += ( ((iy==0)            ? 0 : -rDielectric(ix, iy-1) * leftBoundary.ValueElecFieldX * meshY[iy-1] / 2) // down
						  + ((iy==vertexNumY-1) ? 0 : -rDielectric(ix, iy)   * leftBoundary.ValueElecFieldX * meshY[iy]   / 2) ); // up
		}
	}

	if (rightBoudary.Type == Second)
	{
		ix = vertexNumX - 1;
		for (iy = 0; iy != vertexNumY; ++iy)
		{
			index = positoner.IndexAt(ix, iy);
			rhs[index] += ( ((iy==0)            ? 0 : rDielectric(ix-1, iy-1) * rightBoudary.ValueElecFieldX * meshY[iy-1] / 2) // down
						  + ((iy==vertexNumY-1) ? 0 : rDielectric(ix-1, iy)   * rightBoudary.ValueElecFieldX * meshY[iy]   / 2) ); // up
		}
	}

	if (upBoundary.Type == Second)
	{
		iy = vertexNumY - 1;
		for (ix = 0; ix != vertexNumX; ++ix)
		{
			index = positoner.IndexAt(ix, iy);
			rhs[index] += ( ((ix==0)            ? 0 : rDielectric(ix-1, iy-1) * upBoundary.ValueElecFieldY * meshX[ix-1] / 2) // left
						  + ((ix==vertexNumX-1) ? 0 : rDielectric(ix, iy-1)   * upBoundary.ValueElecFieldY * meshX[ix]   / 2) ) ; // right
		}
	}

	if (downBoundary.Type == Second)
	{
		iy = 0;
		for (ix = 0; ix != vertexNumX; ++ix)
		{
			index = positoner.IndexAt(ix, iy);
			rhs[index] += ( ((ix==0)            ? 0 : -rDielectric(ix-1, iy) * downBoundary.ValueElecFieldY * meshX[ix-1] / 2) // left
				         +  ((ix==vertexNumX-1) ? 0 : -rDielectric(ix, iy)   * downBoundary.ValueElecFieldY * meshX[ix]   / 2) ); // right
		}
	}

	//////////////////// for Dirichlet boundary conditions //////////////////
	if (leftBoundary.Type == First)
	{
		ix = 0;
		for (iy = 0; iy < vertexNumY; ++iy)
		{
			index = positoner.IndexAt(ix, iy);
			rhs[index] = leftBoundary.ValuePotential;
		}
	}

	if (rightBoudary.Type == First)
	{
		ix = vertexNumX - 1;
		for (iy = 0; iy != vertexNumY; ++iy)
		{
			index = positoner.IndexAt(ix, iy);
			rhs[index] = rightBoudary.ValuePotential;
		}
	}

	if (upBoundary.Type == First)
	{
		iy = vertexNumY - 1;
		for (ix = 0; ix != vertexNumX; ++ix)
		{
			index = positoner.IndexAt(ix, iy);
			rhs[index] = upBoundary.ValuePotential;
		}
	}

	if (downBoundary.Type == First)
	{
		iy = 0;
		for (ix = 0; ix != vertexNumX; ++ix)
		{
			index = positoner.IndexAt(ix, iy);
			rhs[index] = downBoundary.ValuePotential;
		}
	}
	//PrintConsole::PrintVector(rhs);
}

void TwoDimProblem::solvePotential()
{
	initRHSVector();
	refreshRHSVector();
	TwoDimSolver twoDimeMath = TwoDimSolver(vertexNumX, vertexNumY);
	twoDimeMath.SolveTwoDim(rhs, potential.DataVector, meshX, meshY, rDielectric, leftBoundary, rightBoudary, upBoundary, downBoundary);
}
